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1 Preliminaries and notation 

Let if he a locally finite point pattern in arising as a realisation of simple point processes 
$ on M'^ [H [9j. In practice, d G {1,2,3}. We shall assume that the points are in general 
quadratic position [H], that is, (a) no d+2 points are located on the boundary of a sphere, and 
(b) in the plane no three points are co-linear; in higher dimensions, no /c+l points lie in a A; — 1 
dimensional affine subspace for k = 2, . . . d. These assumptions are satisfied almost surely 
for realisations of a Poisson process with locally finite intensity function A : M*^ — > [0, oo) or, 
more generally, for Gibbs point processes defined by their probability density with respect to 
such a Poisson process. 

Any point pattern ip gives rise to two interesting tessellations. First consider the set 

C{xi I := {y £ R"^ : \ \xi - y\ \ < \\xj - y\\ ^xj G 93} 

that consists of all points in M°' that are at least as close to Xj G (/? as to any other point of 
If, which is called the Voronoi cell of Xj. The ensemble of all Voronoi cells is the Voronoi 
tessellation of ip [20]. An equivalent definition is 

C{xi \(p) = Pi H{xi,Xj), 

where H{xi,Xj) is the closed halfspace {y G M'^ : (y — {xi + Xj)/2, Xi — Xj) > 0} consisting of 
points that are at least as close to Xi as to Xj. In M^, for Xi < Xj, H{xi,Xj) = {—00, {xi+Xj)/2]. 
In the plane, H{xi,Xj) is the closed halfplane bounded by the bisecting line L{xi,Xj) of the 



segment connecting Xi and xj that contains Xi. Note that the Voronoi cells are closed and 
convex, but not necessarily bounded. 

Under our assumptions, intersections between k = 2, . . . ,d+ 1 different Voronoi cells are 
either empty or of dimension d — k + 1. In particular, 

d+l 

Pi C{xi I (^) ^ <^ b{xi, Xd+i) n 99 = 

i=l 

where b{xi, . . . , Xd+i) is the open ball spanned by xi, . . . Xd+i on its boundary, and in that 
case is a single point, usually referred to as a vertex of the Voronoi diagram. 

Vertices can be used to define the second tessellation of interest to us in this paper, the 
Delaunay tessellation. Indeed, suppose that (p contains at least d+l points. Each Voronoi 
vertex arising as the intersection oi d+l cells C (xj | (p) defines a closed simplex, the convex 
hull of {xi, . . . ,Xrf+i}, which is called a Delaunay cell [5] and denoted by D{xi, . . . ,Xrf+i). 
Note that for d = 1, Delaunay cells are intervals, whilst in the plane they form triangles. 
An alternative, equivalent, edge based construction is to join points xi,X2 € that share 
a common Voronoi border C{xi | if) n C(x2 | 7^ into a Delaunay edge. In this case, 
xi and X2 are called Voronoi neighbours. The set of neighbours of xi in (p is denoted by 
N{xi I (/?). Either way, the partition of space formed by the Delaunay cells is referred to as 
the Delaunay tessellation. The union of Delaunay cells containing Xj G is known as the 
contiguous Voronoi cell W{xi \ ip) of Xj in 93. 




Figure 1: A set of thirty points with their Voronoi (dashed lines) and Delaunay (solid lines) 
tessellations. A contiguous Voronoi cell is indicated by shading. 

For more details, including an historical account, the reader is referred to the compre- 
hensive textbooks |11[ [T2] . An illustration is given in Figure [H which was obtained using the 
DELDIR package [T9]. 
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2 Delaunay tessellation field estimator 



Recently, Schaap and Van de Weygaert [14:\ [T5] proposed to estimate the intensity function 
of a spatial point process by the so-called Delaunay tessellation field estimator (DTFE). 
The method estimates the intensity at points in a realisation reciprocal to the volume of 
their contiguous Voronoi cell, and distributes these estimated field values over Delaunay cells 
by linear (or other) interpolation. They also consider interpolation of fields x i— > f{x) € 
M'^ observed at sampling points. Earlier suggestions to use Voronoi tessellations for field 
interpolation include those by Ord [13] and Sibson [16]. 

Based on extensive simulations, Schaap and Van de Weygaert claim that, in contrast to 
kernel estimators [I], the DTFE preserves the total mass of the field and fine structural details, 
appears to result in smooth interpolation, adapts itself to the local scale and geometry, and 
is relatively robust. The limitations of the method lie in its sensitivity to measurement error, 
boundary effects, and triangular artefacts [15]. Our aim in this paper is a rigorous analysis 
of this estimator. 

Throughout this paper, let $ be a simple point process on M*^ having realisations in general 
quadratic position for which the expected number of points placed in bounded Borel sets is 
finite so that its (first order) moment measure exists as a cj-finite Borel measure. Furthermore, 
assume that the moment measure is absolutely continuous with respect to Lebesgue measure 
with Radon-Nikodym derivative A : M'^ ^ [0, cxd), its intensity function. 

Definition 1. Consider a point process ^ observed in a convex bounded Borel subset A of 
M^. For X e^nA, define 

■= |W^(x%n^)|' 

where \ ■ \ denotes d-volume. For any xq ^ A in the interior of some Delaunay cell, define 

^y-=^i E ^) (2) 

a;G<I>nD{a;o|*nA) 

as the average of the estimated intensity function values at the d+1 vertices x of the Delaunay 
cell D{xq \ ^ f] A) containing xq. 

A few remarks are in order. Should a particular realisation ip of ^ happen to contain 
less than d + 1 points in A, the intensity function estimate may be set to zero, or (cf. the 
Lemma below) to the number of points divided by \A\. On the sides of the Delaunay cells, 
any averaging may be used - it is a null set. Finally, A(a;o) is set to zero for points that do 
not fall in any Delaunay cell. 

Edge effects arise due to the fact that <I> is not observed, only $ n A, the Delaunay 
tessellation of which partitions the convex hull of ^ f] A Q A. Such effects may be dealt 
with in many ways. For example, one might use torus corrections, add arbitrary points on 
the boundary of A (the corners for example in the generic case of a cube), or draw lines 
orthogonal to the edges emanating from points on the boundary of the convex hull, etc. 
Further examples can be found in chapter 6 of |12j . 
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Lemma 1. (Schaap and Van de Weygaert fT4\ 

Let ip be a realisation of the simple point process ^ containing at least d + I points in A. 
Then the estimator of Definition[J\ preserves total mass, that is, 

X{xo)dxo = n{ip n ^4), 

J A 

the number of points of ip in A. 

Proof: Write T>{(p H A) for the family of Delaunay cells defined hy ip D A, and note that 



X{xo) dxo 



Dj(^v{4)nA) 



E 



\W{x \'pnA)\ 



E 

xeipnA 



1 



\w{x I ipnA)\ 



i{xeD,}\D,\ 

Dj(^V{<t>r\A) 



n{ip n A), 



cf. [151 p. 62 ff.]. 



□ 



3 Mean and variance of the Delaunay tessellation field esti- 
mator 

In this section, we derive the first two moments of the Delaunay tessellation field estimator. 
Our first result concerns the expectation. 

Theorem 1. Let ^ be observed in a convex bounded Borel subset A, and, for a point pattern 
ip with n{p f] A) > d + 1 in general quadratic position, set 

5(xo X, (/J := — — — -— , 3) 

\W[x I (/? n 74)1 

for xq ^ A\ip, X ^ ip, and let g{x \ x, p>) := {d + | (^9 n ^)| if x G (j) Ci A. Then the 

Delaunay tessellation field estimator defined by (0j and ([ip has expectation 



E 



A(xo) = / ^xigixQ I x,$)] X{x)dx, 



A 



where E^,. denotes the expectation with respect to the Palm distribution of ^ at x. 

For patterns p with less than d+ \ points falling in A, it is also possible to write X{xq) = 
YlxeipnA 9(^0 I 2;, p>) with the function g chosen to suit the particular type of edge correction 
adopted, see Section [2l 

Proof: Note that _____ 

Kxo) = Y 9{xo I x,$). 
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Hence, by the Campbell-Mecke theorem [18 



EA(xo) 



[g{xo I x,^)] X{x)dx. 



□ 



Recall that the second order factorial moment measure /i^^^ is defined in integral terms 



by 



E 



J2 /(^I'^s) 



f{xi,X2) dfl^'^\xi,X2 



(4) 



for any non-negative measurable function /. The sum is over all pairs of different points. 
We shall say that the second order factorial moment measure exists, if it is locally finite. 
If furthermore absolutely continuous with respect to the 2-fold product measure of 

Lebesgue measure with itself, a Radon-Nikodym derivative exists known as second order 
product density and denoted by p^'^^ . In this case, (HD reduces to 



f{xi,X2) p^'^\xi,X2) dxi dx 



2- 



Theorem 2. Let <I> be observed in a convex bounded Borel subset A and define the function 
g by Assume that second order product densities exist. Then the Delaunay tessellation 
field estimator defined by (0j and (OP has variance 



Var(A(xo)) 



A J A 



Eg [g{xo I X, g{xo \ y, $)] p^^^ (x, y) dx dy 



+ j^E^[g\xo\x,^)]X{x)dx- (^j^E^ 



[g{xo I x,^)] \{x)dx , 



(2) 

where Ex,y denotes the two-fold Palm distribution of^. 



Proof: Remark that 



E 



A(xo) 



E 



^ 9{xo I x,^)g{xo I 

x,y£^nA 



+ E 



.xG^nA 



The cross term on the right hand side is equal to 

^xliaixQ I x,<i>)g{xo I y,$)] p^^\x , y) dx dy , 



AJA 



see e.g. [4J, where Ex,y denotes the two-fold Palm distribution of <I> [8j. Another appeal to 
the Campbell-Mecke theorem yields 



E 



g'^ixQ I x,$) 

.a;G$nA 



/ E, [g\xo\x,^)] AC 
J A 



x) dx. 
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Finally, the variance is obtained using Theorem [TJ 



□ 



In general, the integrals involved in Theorems [TH2] must be evaluated by numerical or 
simulation methods. 



4 Comparison to a classic estimator 

The classic estimator of intensity is the kernel estimator 

.-—r\ n{<^nb{xo,h)nA) 
\b{xo,h)nA\ ' 

proposed by Berman and Diggle [Ij. The estimator can be regarded as a kernel estimator |17j 
with kh{xQ I x) = l{||x — xoll < /i}/|6(xo, /i) n^l, where b{xo,h) denotes the open ball around 
xq with radius h > 0. The choice of bandwidth h determines the amount of smoothing. 

Note that when the bounded observation window A 7^ is open, one never divides by 
zero. In fact, a stronger statement can be made. The function x 1— > \b{x, h) nA\ is continuous 
and attains its minimum on the closure A. Since any point on the boundary dA always has 
a neighbour within distance h in A, inix^A \b{x, h) f] A\ > 0. Further details may be found 
e.g. in m EKE]. 

Although ([5]) is a natural estimator, it does not necessarily preserve the total mass in A 
[15] . nor is it based on a generalised weight function [T7]. It is not hard to modify the edge 
correction in ^ to define an estimator [10] that does preserve total mass and is based on a 
weight function. 

Definition 2. Consider a point process $ observed in an open bounded Borel subset A ofW^ . 
For xq S a, define 

^^^"°)^= \b{xh)^A\ • 

Lemma 2. The estimator of Definition\^ is a generalised weight function estimator with 
kernel kh{xo \ x) = l{\\x — xo|| < h}/\b{x^h) n A\ that preserves total mass, that is, 



A 



XKixo) dxo = n A), 



the number of points of ^ in A. 
Proof: Note that 

kh{xo I x) dxo = [ — dxo = 1 
J A \b{x,h)nA\ 

for all X A, that is, Ax(') is a generalised weight function estimator. Furthermore, for any 
realised point pattern 93, the restriction (pCi A in A is finite and 



E 



l{||x-xol| < h} 



\b(xh)nA\ 



d^o= ^ A \b{xh)nA\ dxo = n{^nA). 



xG^pnA 



A 
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□ 

Note that the Delaunay tessellation field estimator is based on an adaptive kernel Q as 
it depends on the underlying point pattern. Indeed, for every x (z A, 



El 



i{xoeD^.-xeD,} 



g{xo X, 4>) dxQ = / — — — dxo = 1. 

A J A \W{x\ipnA)\ 

A clear advantage is that the problem of choosing the bandwidth is avoided. 

In order to assess the quality of the estimator, we proceed to compute its mean and 
variance. 

Theorem 3. Let <I> be observed in a bounded open Borel subset A. Then, the estimator of 
Definition [H has expectation 



E 



Ak(xo) 



l{x G b{xQ, h)} 



I A \b[x,h)r\A\ 
Proof: By the Campbell-Mecke theorem 



A(x) dx. 



E 



E 



\b{xh) r\A\ 



l{||3:-a:o|| < h} 
\b{x,h) r\A\ 



\{x) dx. 



□ 



If we compare Theorem [3] to Theorem [H the Palm expectation E^. [g{xQ j x, ^] is replaced 
by kh{xQ I x), as the latter does not depend on the point process <I>. 

Theorem 4. let <I> be observed in a bounded open Borel subset A and assume that second 
order product densities exist. Then 



Var(Ax(xo)) 



p(2)(x,y)-A(x)A(y) 



'(fe{xo,h)nA)2 \b{x,h)r\A\ \b{y,h)nA\ 
Proof: Regarding the second moment, note that 



dx dy + 



/ 

J b(xo 



A(x) 



h)nA \b{x,h)nA\' 



dx. 



E 



= 4 E 

\x,y£^nA 



l{||x-xo|| < h} l{\\y-XQ\\ < h} 

\b{x,h)nA\ \b{y,h)nA\ 

l{l|x — xoll < /i} 

|6(x,/i)nA|2 



Then rewrite the expectations as integrals with respect to p*-^^ and A respectively to obtain 
that the variance of Ai^(xo) is equal to 



/ / 

Jb{xo,h)nA Jb(xo 



h)nA \b{x,h)nA\ \b{y,h)nA\ 



p^'^\x, y) dxdy + / 

Jb{xo 



A(x) 



b(xo,h)nA |6(x,/i) n A|2 



dx. 
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An appeal to Theorem [3] completes the proof. 



□ 



The result should be compared to that of Theorem [21 

Similar arguments as those in the proofs of Theorems [3] and U] applied to the classic 
Berman-Diggle estimator ^ give mean 



1 



\b{xo,h) nA\ 



X{x) dx 



b{xo,h)nA 



and variance 



1 



\b{xo,h) nA\' 



/ A(x) dx + 

lb{x(),h)nA) J(b{xo,h)nA)'2 



j 

J(b\ 



p^'\x,y)-X{x)X{y) 



dx dy 



Note that for xq AQ 6(0, 2h) separated by 2h from the boundary of A, no edge correction 
is necessary, and both kernel estimators are identical. 

The disadvantage of kernel estimators is that they involve a bandwidth parameter h; the 
larger h, the smoother the estimated intensity function. For specific models, h may be chosen 
by optimisation of the (integrated) mean squared error [6]. In practice, in a planar setting, 
Diggle [6] recommends to choose h proportional to n~^/^, where n is the observed number of 
points. For a fixed bandwidth, neither the Berman-Diggle estimator nor the modification of 
Definition [2] is universally better. For examples, the reader is referred to [10\ - 



5 Intensity estimation for Poisson point processes 

In general, the integrals involved in Theorems [iHll have to be evaluated numerically. An 
exception is the case where $ is a Poisson point process with a locally finite intensity function. 

Corollary 1. Let ^ be a Poisson point process observed in a convex bounded Borel subset A. 
Then, 



E 



E [g{xQ I X, <I> U {x})] \{x) dx 



and 



Var(A(xo)) 




E [g{xQ I X, ^> U {x, y}) g{xQ | y, ^> U {x, y})] \{x) \{y) dx dy 



A J A 

E [g^{xQ I X, $ U {x})] A(x) dx 



E [g{xQ I X, $ U {x})] A(x) dx 



Proof: For a Poisson process, the Palm distribution at x is equal to the superposition of its 

(2) 

distribution P with an extra point at x, the two-fold Palm distribution Pi^^ is the superpo- 
sition of P with X and y. Furthermore, p^'^\x,y) = A(x) A(y) is a product density. Plugging 
these results into the expressions of Theorems [THl] completes the proof. □ 
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Corollary 2. let ^ be a Poisson point process observed in a bounded open Borel subset A 
and assume that second order product densities exist. Then, 



Var(Ai^(xo)) 



A(x) 



lb{xo,h)nA \bix,h)nA\'^ 
Proof: Use that p^'^\x,y) = X{x) X{y) and apply Theorem 21 



dx. 



□ 



The variance of the Berman-Diggle estimator is J^j-^^ ^^^^^ X{x) dx/\b{xo, h) f] . 

For stationary Poisson processes, even more can be said. In the remainder of this section, 
define g as in ^ with A = M"'. 

Theorem 5. Let ^ be a stationary Poisson point process in with intensity A > 0. Then, 
the Delaunay tessellation field estimator A(0) is asymptotically unbiased. 

Proof: Let yi, . . . , y^) be the open ball spanned by the points x, yi,...,yd on its 
topological boundary, and let D°{x, yi, . . . , yd) be the open simplex that is the interior of the 
convex hull of {x, yi, . . . , y^}. Recall that the points x, yi, . . . , y^ define a Voronoi vertex, or, 
equivalently, a Delaunay cell if and only if there are no points in b{x, yi, . . . , ya). 
By Corollary [H asymptotically 



E 



A(0) 



A / E[y(0 [ x,$U{x})] dx 



A / E 



A / E 



A / E 



A / E 



E 

{yi,---,yd}c'S> 

E 

{yi,---,yd}c<s> 

E 

{zi,...,za}C<S>-^ 

E 

{zx,...,Zi\<Z^ 



1{0 £ D°{x, yi, ■ ■ ■ , yrf); 6(x, yi, . . . , y^) n U {x}) = 0} 
|l^(x I «>U{x})| 

1{0 elj°(x,yi,...,yrf);6(x,yi,...,yrf)n^ = 0} 
|T^(x I $U{x})| 

l{-x £ i^°(0, zi, . . . , zs)\ 6(0, zi, . . . , Zrf) n = 0} 

1^(0 I $-..u{o})l 

l{-x e D°(0, zi, . . . , Zd); 6(0, zi, . . . , Zrf) n $ = 0} 



dx 



dx 



|VF(0 I ^U{0})I 



dx 



by stationarity. Hence, by Fubini's theorem, 
E A(0) = AE ' ^ 



Et,...,.,}c* l^°(0' ^1, • • • > 1^0, zi, . . . , Zd) n cl> = 0} 



= AE 



|PF(0 I «>U{0})| 
|M/(0 I ^U{0})[ 



|Ty(0 I $U{0})| 



A. 



□ 



9 



The asymptotic variance of the Delaunay tessellation field estimator increases quadrati- 
cally with A with a constant multiplier that depends on the dimension. The proof rests on 
the following two lemmata. 

Lemma 3. Let ^ be a stationary Poisson point process in M.'^ with intensity A > 0. Then, 



C{X,d) 



E [g{0 I X, $ U {x, y}) g{0 [ y, $ U {x, y})] A(x) A(y) dx dy 

\w{<d I $u{o,x}) npr(x I «>u{o,x})r 



E 



dx. 



|VF(0 I $U{0,j;})| \W{x I «>U{0,x})| 
where Ei denotes expectation with respect to a unit intensity Poisson point process. 
By the Nguyen-Zessin formula [9j, alternatively 



C{\d) 



AE 



|VF(0 I $U{0})| 



E 



S/GA^{0|<I>U{0}) 



|VF(o I $u{o}) nTy(y [ $u{o})| 

\W{y\^\j{m 



A^ El 



\W{Q I $U{0})| 



E 



|VF(0 I $U{0}) nT^(y I $U{0})[ 



yeM{o\^vj{Ci}) 



\W{y I <I>U{0})| 



Proof: Write ^d-i for sets of d — 1 distinct points in Then, as A(x) = A is constant, 
and 5^(0 | x, <I> U {x, y}) g{0 | y, $ U {x, y}) vanishes when x and y do not belong to the same 
Delaunay cell containing in its interior. 



C{X,d) = A^ 



E 



1{0 € D°{x, y, z); b{x, y, z) n ^ = 0} 
^ \W{x\^U{x,y})\\Wiy \<i>U{x,y})\ 



dx dy 



AM / E I 

Because of stationarity. 



E 



l{-x E D°{0, y -X, z); 6(0, y - x, fl 



0} 



\W{0 I U {0, y - x})| \W{y - x | «>_^ U {0, y - x})| 



dx dy. 



C{X,d) 



E 



E 



E 



l{-x e D°{0, y - X, z);b{0, y - x, z) n $ = 0} 
\W{0 I ^>U{0,y-x})j |M^(y-x [ $U{0,y-x})| 



dx dy 



dx dy. 



|T^(0|$U{0,y})||VF(y|<I>U{0,y})| 
Scaling by X^^"^ yields that A~^ C{X,d) is equal to 

'Eze<fd-i H-^^^'^x e D°{0, X^'^y, X'^/'^z); 6(0, A^^^y, X^/'^z) n A^/^^* = 0} 



E 



A-i|VF(0 I Ai/««$u{0,AV'^y})|A"i|VF(AV'^y | AV'i^ u {0, Ai/<iy})| 



dx dy. 
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bmce is a unit intensity Poisson point process, we obtain 

"E.g$,_i e ^°(o, y, zy,b{o, y,z)n<^ = 0} 



A-^C(A,(i) 



El 



\W{0 I cDu{0,y})||T^(y | cDu{0,y})| 
An appeal to Fubini's theorem to integrate out over x completes the proof. 



dx dy. 



□ 



Lemma 4. Let ^ be a stationary Poisson point process in M'^ with intensity A > 0. Then, 



C'{X, d):= E [^2(0 I X, $ U {x})] X{x) dx = A^ Ei 



1 



\W{0 I «>U{0})[ 

where Ei denotes expectation with respect to a unit intensity Poisson point process. 
Proof: Using A(x) = A and argueing as in the proof of Theorem [5l we get 

7^ 



C'{\,d) = A / E 



A / E 



E 

{zi,...,za}C<S> 

E 

{zi,...,za}C'l> 



l{-x e D°iO, zi,..., zd);b{0, zi, . . . , z^) n ^ = } 
\W{0 I $U{0})| 

l{-x G D°{0, zi,..., Zd); b{0, zi, . . . , z^) n $ = 0} 



dx 



\W{0 I $U {0})|2 



(7) 



as —X belongs to a single Delaunay interior. Write for sets of d distinct points in $ and 
scale each point in ^ by X^^^ to obtain that C'{X, d) is equal to 



A / E 



E 



l{-X^/'^x G D°{0, aV'^z); 6(0, AiZ-^z) n AV'^^ = 0} 
A-2|Ty(0 I Ai/rf$U{0})|2 



dx 



which, since is a unit rate Poisson process reduces to 

"E{.i,...,^,}c* G 2?°(0, zi, . . . , Zrf); 6(0, zi, . . . , Zd) n $ = 0} 



A^ / El 



\W{0 I «5U{0})[ 



dx. 



The sum of d- volumes of Delaunay cells involving is that of its contiguous Voronoi cell, and 
we conclude that 



C'{X,d) = A^Ei 



1 



\W{0 I $U{0})| 



□ 



The above results can be summarised as follows. 

Theorem 6. Let ^ be a stationary Poisson point process in with intensity A > 0. Then, 
the Delaunay tessellation field estimator A(0) has asymptotic variance c^X^ with 

1 L |VF(0 [ $U{0})nT^(y I $U{0})[ 



Cd = Ei 



\W{0 I $U{0})| 



1 + 



E 



j/eA^(o|*u{o}) 



\Wiy I $U{0})| 



1. 
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Note that the classic Berman-Diggle estimator ([5]) is asymptotically unbiased with vari- 
ance Ao;^^ h~'^, where cod is the volume of the unit ball in M*^. In words, the Berman-Diggle 
estimator is more efficient whenever the average number of points per test ball exceeds I/q. 



6 Poisson processes on the line 

For one-dimensional Poisson processes, the distribution of the contiguous Voronoi cell can be 
calculated explicitly for arbitrary intensity functions. For simplicity, assume that A = [—w, w] 
is an interval of radius w > either side of the origin. 
The following lemma is well-known. 

Lemma 5. Let ^ be a Poisson point process on [—w,w] with finite intensity function A and 
write A{a,b) = jj^ X{x) dx for the moment measure of {a,b) for any —w < a < b < w. For 
X G {—w,w), define the random variables 

^^{x) := max{y e {-w}U{^n[-w,x))}; 
<i)+(x) := mm{y e {w} U n {x,w])}. 

Then, their distribution functions are given by 

F~{t) = exp [-A{t,x)] 

fort G {—w,x), with an atom of mass P{^^{x) = —w) = exp [—A{—w,x)] at —w, respectively 

F+{s) = 1 - exp [-A(a;,s)] 

for s € {x,w) with an atom at w of mass P{^^{x) = w) = exp [— A(x, u;)] . Moreover, for 
fixed X, <I>^(x) and ^~{x) are independent random variables. 



6.1 Expectation of the DTFE for Poisson processes on the hne 

Note that on the real line, the contiguous Voronoi cell W{x \ {^L){x})r][—w, w]) is the interval 
Thus, Lemma [5] can be used to calculate the moments of the Delaunay 
tessellation field estimator. In this section, we shall deal with edge effects by placing two 
ghost points at the borders —w and w. 

Theorem 7. Let ^ be a Poisson point process observed in A = [—w, w] for some w > with 
locally finite intensity function A : M ^ [0, oo). Then, for xq € A, 



E 



X{xo) 



XQ 



W J Xq 



A{t,s)\{s)\{t) ^_^(,^,) 
s - t 



dt ds + 



A{-w,w)e~^'^-'"''"^ 



2w 



Ai-w,s)X{s) ^_A(„^,,) ^ r« A{t,w)Xit) ^_A(,,^) 



XQ 



W + S 



dt. 



w 



(8) 
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Proof: Fix xq 7^ x G {—w,w), and let (f he a realisation of which we augment by —w 
and w in order to obtain bounded Delaunay cells. Since almost surely, x ^ $ and xq ^ 
assume xq, x ^ ip, and consider g{xQ \ x,if[J {x}) as defined in ([3]). Note that xq belongs to a 
single Delaunay cell interior. If x is no endpoint of this cell, g{xo \ x, ipL){x}) = 0. Otherwise, 
g{xQ \ x,(p[J {x}) = !/((/? ^"(xo) — (p~{xq)), cf. Lemma O 

First, assume x < xq. By Lemma [5] applied to the point x, 



E [g{xo I X, $ U {x})] 



dF-(t)dF+(s) e-'^(-"''"') 



w J XQ s t 2w 



+ r iifL e-'^(-'^) ds + r ^ e-^(*'-) dt + r r ^^^^ e-^^*-^) dt ds. 
Jxo w + s 

Similarly, for xq < x. 



XQ W -\- S J ~W ^ ^ J —W J Xq ^ ^ 



w + s 7-1/1 w — t j _yj s — t 



dt ds. 



By Theorem [H the expectation of the Delaunay tessellation field estimator is as stated for 
Xo € {—w, w). 

It remains to consider xq = —w or w. In the first case, ip~{xo) is replaced by —w; for 
Xo = w, ip~^{xq) becomes w in the evaluation of g{xo \ x,ipU {x}). Thus, for example, 

E I $ u = e-«-«> r = + r d,, 

Jx w + s 2w Jx w + s 

with a similar expression for xo = w. Upon integration, ([8j) is obtained, under the convention 
that integrals over intervals of zero length vanish. □ 

In general, ([8]) must be evaluated numerically. For the homogeneous Poisson process, 
analytic evaluation is possible. In fact, it can be shown that the estimator is unbiased even 
near the borders of the observation interval. 

Corollary 3. Let ^ be a stationary Poisson point process observed in A = [—w, w] for some 

w > with intensity A > 0. Then, the Delaunay tessellation field estimator A(xo) is unbiased 
for all xq G A. 

Proof: For a stationary Poisson point process, the double integral in dHj) reduces to 

and in particular vanishes for xq = —w or w. The three border correction terms are equal 
to Ae-^^"", to Ae-^"'(e-^^o - e"^""), and to Ae-^"'(e-^^o - e"^"'), respectively. The sum of all 
four terms is A, so the estimator is unbiased. □ 

Note that the Berman-Diggle estimator is unbiased as well, but that this may not be true 
for dnj) due to edge correction near the border. 
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6.2 Variance of the DTFE for Poisson processes on the hne 

In this section, we derive the asymptotic variance of the Delaunay tessellation field estimator 
for a stationary Poisson process on the line. The result can be used to approximate the 
variance when the underlying intensity function is smoothly varying. 

Theorem 8. Let ^ be a stationary Poisson point process observed in A = [—w, w] for some 
w > with intensity A > 0. Then, as w ^ oo, the Delaunay tessellation field estimator A(0) 
has asymptotic variance 2A^(2 — vr^/G) ~ 0.7A^. 

The result should be compared to \/{2h) for the Berman-Diggle kernel estimator [1], see 
also [lOj. If 2Xh > 1.4, that is the average number of points per bin at least 1.4, kernel 
estimation is the better choice. Naturally, in order to compute A(xo), two points of the 
underlying process are used. 

In order to give the proof, some special function theory is needed. Let x > 0. Recall that 
the exponential integral is defined as 

/•CO ^—tx t-OO g— 11 

E,(x)=y^ ^dt = j^ —du. 

Its integral satisfies 

E2{x) = / Ei{s)ds = e"^ - xEi{x). 

J X 

In the limit, -E'i(O) = oo and £'2(0) = 1. Furthermore, 

2 



foo 2 

/ ue"" Ei{uf du = 2 

Jo 6 



See for example ^IJ for further details. We shall also need the equation 

7 + log(ac) + e°-'^ El (ac) 



Te"^ Ei{ax)dx 
Jo 



a 

where a and c are strictly positive constants, and 7 ~ 0.577 is the Euler-Mascheroni constant. 

Proof: By Theorem [3l asymptotically E A(0) = A. For the variance, by Theorem [21 we 
need to evaluate two further integrals. Now, argueing as in the proof of Theorem [TJ 

/ E[g\xo\x,'!>U{x})]X{x)dx= H T AM^if)^ e-A(M 

J A J -w Jxo ^) 

A{-w,w)e-^(-^-) r A{-w,s)X{s) ^^(_,,) n A{t,w)X{t) , (.^^ 
Since the intensity function is constant and we took xq = 0, ([9|) reduces to 

\„— 2Aui fw \„— As fO \„Xt rO fw \3„At_— As 

AC , / Ae \-Xw I Ae ,. / / A e e 



+ Xe~^'" / ds + Xe-^"" / dt + / / dtds. 

2w Jo w + s J_^ w - 1 Jq s-t 



14 



Clearly, the first term above converges to as ^ oo. Due to symmetry, the two middle 
terms are equal. Note that 

f-w \ —\{s+w) r2\w p — u 

2A / ds = 2X^ du = 2X^ [EUXw) - Ei(2Xw)] , 

Jo s + w Jxw u 

which converges to zero as w ^ oo. Moreover, 



-oo JO 



dtds = AM Ei{-Xt) dt = X^E2{<d) = X^ 



To calculate the double integral in Theorem [2l let x 7^ y be points of {—w^w)., fix xq 
{x, y, — tt;, It;}, and let 99 be a realisation of which we augment by —w and w in order 
to obtain bounded Delaunay cells. Since almost surely none of x, y or xq lie in assume 
xo,x,y ip, and consider g{xo \ x,ifU {x,y}) as defined in ([3]). Note that xq belongs 
to a single Delaunay cell interior. If x and y are not both endpoints of this cell, g{xo \ 
x,ifU {x, y}) g{xo \ y,ipU {x, y}) = 0. Otherwise, without loss of generality, x < xq < y, and 
y(xo I x,ipU{x,y}) = l/(y - 92" (xq)) and y(xo | y,</? U {x,y}) = !/((/?+ (xq) - x). 

Thus, for X < Xq and y > xq, let F~ and be the cumulative distribution functions of 
$~(xo) and <^'*'(xo). By Lemma El 



E [y(xo I X, $ U {x, y}) y(xo | y, $ U {x, y})] 



dF-{t)dF+{s) 
y (y -t)is- x) 



X{s) 



y iw + y){s- x) 

p-A{-w,w) 



-A{-w,s) 



+ 



m 



(y -t){w - x) 



-A(i,«,) 



+ 



By symmetry, 



{w + y){w- x) J_^ Jy (y -t){s- x) 
E [y(xo I X, $ U {x, y}) y(xo | y, ^> U {x, y})] A(x) A(y) dx dy 



AJA 



+ 2 
+ 2 
+ 2 



rw 

J XO 



X{t)e 



w — X 
-A{-w,s) 



dy 

XO w + y 

^0 A(x) , 

dx 

s — X 



W) '-' J Xq 

XO x(x) , 
dx 



w — X 



w + y 



dy 



ds 



y-^ 



dy 



A(t)A(s)e-^(*'^) 



A(x) 



dx 



s — X 



dt 

Ky) 

^oV-t 



dy 



dt ds. 



(10) 



For Xq G {—WjW}, formula (jlOp holds true under the convention that integrals over intervals 
of zero length vanish, as in this case xq cannot belong to any Delaunay cell with endpoints 
X < xo < y. 
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Next, we plug in xq = and A(-) = A, and consider each integral in (jlOp in turn. The 
main term is the four fold integral 



/■■(« fX pw gAtg— As 



w JO J —w J y 



{y-t) {s-x) 



dx dy dt ds. 



Its limit as w ^ oo is 

_r>0 /-oo 

2A^ 



-oo JO 



■ dt / ds 



y-t 



s — X 



dx dy 



(u) dx du 



/U fco pv poo 

/ e^^y-'''^Ei{\{y-x)fdxdy = 2X^ / e"Si(i 
ooJo J —oo J —Xx 

2A2 / / e''Ei{ufdydu = 2\^ j u e"" Ei{uf du = 2 X"^ {2 - tt'^ / 6) , 

Jo Jy Jq 

upon a change of integration order. 

The first term in (jlOh reduces to 2e^^^'"(A log 2)^ for a homogeneous Poisson process, 
which tends to zero as t(; — > oo. 

It remains to consider the sum of the two three fold integrals in (|10p 



rw f-w g— A(s-l-io) 



Jy is-x){y + w) 



dx dy ds 



which can be written as 

dy \ e-^(^+'") 



4A^ 



w rw / rs 



Jo \Jo y + wj s + x 



dxds <4 A^ log 2 / e 



w g-A(s+x-) 
S + X 



ds dx 



4 AMog 2 / e 



-Xw+Xx 



El {Xx) - El (Ax + Xw)] dx = 4:X'^h{X,w) log 2, 







where 



/i(A, w) 



-Xw 



/•Xw 

/ e"[-Ei(u) --Ei(n + Aw)](iu 
Jo 



/■Au) i-ZAw 

/ e" ^1 (n) du - e-^^"" / Ei (n) 

= €-^""-1 + (e"^'" + e-^^"') log(Au;) - e'"^^"" \og{2Xw) 
+ ^i(Au')(l + e-^"')-^i(2Au') 



(■2Xw 



tends to zero as t/; — > oo. The proof is finished upon collection of all terms. 



□ 



As a corollary, the proof gives an expression for the second moment of the Delaunay 
tessellation field estimator of the intensity function for Poisson processes with not necessarily 
constant locally finite intensity function on intervals of the form [—w^ w\ by combining ([9])- 
([10]). A slightly simpler proof can be obtained by an appeal to Theorem O but such a proof 
cannot be generalised to non-homogeneous Poisson processes. 
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7 Discussion 



In this paper, we analysed Schaap and Van de Weygaert's Delaunay tessellation field estimator 
[14:\ [TS] for the intensity function of a point process. We expressed its mean and variance 
in terms of the first and second order factorial moment measures of the underlying point 
process, and placed the estimator in the context of adaptive kernel estimation. We then 
focussed on Poisson point processes, and showed that for stationary Poisson processes, the 
DTFE is asymptotically unbiased with a variance that is proportional to the squared intensity. 
The proportionality constant depends on the dimension. For d = 1, explicit calculation is 
possible. For d = 2, we used the DELDIR package [T9] to obtain C(A,2) « 0.8A^ and 
C"(A,2) ~ 0.6A^, see Lemma [3] and [H Note that in the plane it is possible to write mean 
and variance as repeated integrals in the spirit of Calka [2], but explicit evaluation seems 
difficult. Simulations for the case d = 3 of most interest to cosmologists can be found in 
Schaap's Ph.D. thesis [15|. 
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